Making maps in R:

Using the sf and tmap Packages

Suggested Citation:

Mendez C. (2020). Making maps in R: Using the sf and tmap Packages. R Studio/RPubs. Available at https://rpubs.com/quarcs-lab/tutorial-maps-in-r

This work is licensed under the Creative Commons Attribution-Share Alike 4.0 International License.

1 Libraries

knitr::opts_chunk$set(echo = TRUE, warning=FALSE)

library(tidyverse)  # Modern data science workflow
library(sf)         # Simple features for R
library(tmap)       # Thematic Maps
library(tmaptools)  # Thematic Maps Tools
library(RColorBrewer) # ColorBrewer Palettes
library(leaflet)    # Interactive web maps
library(rgdal)      # Bindings for the Geospatial Data Abstraction Library
library(rgeos)      # Interface to Geometry Engine - Open Source 


# Change the presentation of decimal numbers to 4 and avoid scientific notation
options(prompt="R> ", digits=4, scipen=999)

2 Tutorial objectives

• Load spatial data files into R
• Join non-spaital data to spatial data files
• Create simple choropleth maps

3 Orginal data sources

The non-spatial datafile is from:

• Öcal, N., & Yildirim, J. (2010). Regional effects of terrorism on economic growth in Turkey: A geographically weighted regression approach. Journal of Peace Research, 47(4), 477-489.

The spatial (shapefile) is from:

• Matthew C. Ingram

4 Import the data

4.1 Non-spatial data

dat <-read.csv("DATASET.csv")

4.1.1 Explore the data

glimpse(dat) 

Observations: 67
Variables: 7
$ id      <int> 1, 2, 3, 4, 5, 7, 8, 9, 10, 11, 12…
$ growth  <dbl> 0.0254, -0.0254, 0.0773, 0.2583, 0…
$ base    <dbl> 6.157, 5.900, 5.896, 5.315, 5.923,…
$ T       <dbl> 3.585, 2.896, 2.763, 3.140, 2.785,…
$ E       <dbl> 5.216, 3.458, 4.995, 3.478, 4.805,…
$ G       <dbl> 6.192, 6.048, 5.950, 5.354, 5.940,…
$ Coastal <int> 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0…

4.1.1.1 Definitions of variables

• growth: growth rate of per capita real income 1987-2001
• base: lograrithm of the per apita income in the base year 1987
• T: average terrorism index
• E: average years of schooling
• G: real per capita government expenditures in 1987
• Coastal: dummy variable which takes the value of one if the province is a coastal province.

4.2 Spatial data

mapData <- read_sf("MAP.shp")

4.2.1 Explore the data

glimpse(mapData)

Observations: 67
Variables: 21
$ ObjectID   <int> 1946, 1949, 1952, 1957, 2116, 2…
$ NAME       <chr> "Adana", "Adiyaman", "Afyon", "…
$ COUNTRY    <chr> "Turkey", "Turkey", "Turkey", "…
$ ISO_CODE   <chr> "TR01", "TR02", "TR03", "TR04",…
$ ISO_CC     <chr> "TR", "TR", "TR", "TR", "TR", "…
$ ISO_SUB    <chr> "01", "02", "03", "04", "05", "…
$ ADMINTYPE  <chr> "Province", "Province", "Provin…
$ DISPUTED   <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
$ NOTES      <chr> NA, NA, NA, NA, NA, NA, NA, NA,…
$ AUTONOMOUS <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
$ COUNTRYAFF <chr> "Turkey", "Turkey", "Turkey", "…
$ CONTINENT  <chr> "Asia", "Asia", "Asia", "Asia",…
$ Land_Type  <chr> "Primary land", "Primary land",…
$ Land_Rank  <int> 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5…
$ Shape_Leng <dbl> 8.589, 5.355, 6.698, 6.553, 4.8…
$ Shape_Area <dbl> 1.6029, 0.8622, 1.4795, 1.1537,…
$ idnum      <int> 1, 2, 3, 4, 5, 7, 8, 9, 10, 11,…
$ province   <chr> "Adana", "Adiyaman", "Afyon", "…
$ id         <dbl> 1, 2, 3, 4, 5, 7, 8, 9, 10, 11,…
$ name_esri  <chr> NA, NA, NA, NA, NA, NA, NA, NA,…
$ geometry   <POLYGON [°]> POLYGON ((36.44 38.22, …

• Check the Coordinate Reference System

st_crs(mapData)

Coordinate Reference System:
  EPSG: 4326 
  proj4string: "+proj=longlat +datum=WGS84 +no_defs"

5 Transform the data

No need to transform any data because both datasets share a common variable id

6 Merge the data

dat_map <- inner_join(
  dat,
  mapData,
  by = "id"
)

6.1 Keep data as sf object

• Keep the data as sf class, so we will not lose the coodinate system

dat_map <- st_as_sf(dat_map)
st_crs(mapData)

Coordinate Reference System:
  EPSG: 4326 
  proj4string: "+proj=longlat +datum=WGS84 +no_defs"

glimpse(dat_map)

Observations: 67
Variables: 27
$ id         <dbl> 1, 2, 3, 4, 5, 7, 8, 9, 10, 11,…
$ growth     <dbl> 0.0254, -0.0254, 0.0773, 0.2583…
$ base       <dbl> 6.157, 5.900, 5.896, 5.315, 5.9…
$ T          <dbl> 3.585, 2.896, 2.763, 3.140, 2.7…
$ E          <dbl> 5.216, 3.458, 4.995, 3.478, 4.8…
$ G          <dbl> 6.192, 6.048, 5.950, 5.354, 5.9…
$ Coastal    <int> 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0…
$ ObjectID   <int> 1946, 1949, 1952, 1957, 2116, 2…
$ NAME       <chr> "Adana", "Adiyaman", "Afyon", "…
$ COUNTRY    <chr> "Turkey", "Turkey", "Turkey", "…
$ ISO_CODE   <chr> "TR01", "TR02", "TR03", "TR04",…
$ ISO_CC     <chr> "TR", "TR", "TR", "TR", "TR", "…
$ ISO_SUB    <chr> "01", "02", "03", "04", "05", "…
$ ADMINTYPE  <chr> "Province", "Province", "Provin…
$ DISPUTED   <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
$ NOTES      <chr> NA, NA, NA, NA, NA, NA, NA, NA,…
$ AUTONOMOUS <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
$ COUNTRYAFF <chr> "Turkey", "Turkey", "Turkey", "…
$ CONTINENT  <chr> "Asia", "Asia", "Asia", "Asia",…
$ Land_Type  <chr> "Primary land", "Primary land",…
$ Land_Rank  <int> 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5…
$ Shape_Leng <dbl> 8.589, 5.355, 6.698, 6.553, 4.8…
$ Shape_Area <dbl> 1.6029, 0.8622, 1.4795, 1.1537,…
$ idnum      <int> 1, 2, 3, 4, 5, 7, 8, 9, 10, 11,…
$ province   <chr> "Adana", "Adiyaman", "Afyon", "…
$ name_esri  <chr> NA, NA, NA, NA, NA, NA, NA, NA,…
$ geometry   <POLYGON [°]> POLYGON ((36.44 38.22, …

7 Plot Thematic Maps

7.1 Quick Map

• A quick map of the Terrorism variable

qtm(dat_map, fill = "T")

The shape dat_map is invalid. See sf::st_is_valid

7.2 Simple Map

tm_shape(dat_map) + tm_fill("T") 

The shape dat_map is invalid. See sf::st_is_valid

7.3 Remove frame

tm_shape(dat_map) + 
  tm_fill("T") +
  tm_layout(frame = FALSE)

The shape dat_map is invalid. See sf::st_is_valid

7.4 Set legend outside

tm_shape(dat_map) + tm_fill("T") +
  tm_layout(legend.outside = TRUE, frame = FALSE)

The shape dat_map is invalid. See sf::st_is_valid

7.5 Set legend size

tm_shape(dat_map) + tm_fill("T") +
  tm_layout(legend.outside = TRUE, frame = FALSE) +
  tm_legend(text.size = 0.3) 

The shape dat_map is invalid. See sf::st_is_valid

7.6 Set color palette

See color reference here

tm_shape(dat_map) + 
  tm_fill("T", palette = "Greens") +
  tm_layout(legend.outside = TRUE, frame = FALSE)

The shape dat_map is invalid. See sf::st_is_valid

tm_shape(dat_map) + 
  tm_fill("T", palette = "-Greens") +
  tm_layout(legend.outside = TRUE, frame = FALSE)

The shape dat_map is invalid. See sf::st_is_valid

tm_shape(dat_map) + 
  tm_fill("T", palette = "viridis") +
  tm_layout(legend.outside = TRUE, frame = FALSE)

The shape dat_map is invalid. See sf::st_is_valid

7.7 Add polygon names

tm_shape(dat_map) + 
  tm_fill("T", palette = "viridis") +
  tm_layout(legend.outside = TRUE, frame = FALSE) +
  tm_text("province", size = "Shape_Area",  auto.placement = F, legend.size.show = FALSE) 

The shape dat_map is invalid. See sf::st_is_valid

Documentation

7.8 Add scale bar

tm_shape(dat_map) + 
  tm_fill("T", palette = "viridis") +
  tm_layout(legend.outside = TRUE, frame = FALSE) +
  tm_text("province", size = "Shape_Area",  auto.placement = F, legend.size.show = FALSE)  + 
  tm_scale_bar(position = c("RIGHT", "BOTTOM"))

The shape dat_map is invalid. See sf::st_is_valid

7.9 Set color intervals

Enter “style =” followed by one of the options below.

• equal: divides the range of the variable into n parts.
• pretty: chooses a number of breaks to fit a sequence of equality spaced ‘round’ values.
• quantile: equal number of cases in each group
• jenks: looks for natural breaks in the data
• cat: when the variable is categorical

tm_shape(dat_map) + 
  tm_fill("T",
          style = "quantile",
          palette = "Reds"
    ) +
  tm_layout(
    legend.outside = TRUE,
    frame = FALSE)

Change the number of intervals in the color scheme and how the intervals are spaced. Changing the number of intervals n = 7. So, we have 7 shades instead of the default 5.

tm_shape(dat_map) + 
  tm_fill("T",
          style = "quantile",
          n = 7,
          palette = "Reds"
    ) +
  tm_layout(
    legend.outside = TRUE,
    frame = FALSE)

7.10 Add histogram

tm_shape(dat_map) + 
  tm_fill("T",
          style = "quantile",
          n = 5,
          palette = "YlOrBr",
          legend.hist = TRUE
    ) +
  tm_layout(
    legend.outside = TRUE,
    frame = FALSE)

7.11 Add borders

You can edit the borders of the shapefile with the tm_borders() function which has many arguments. alpha denotes the level of transparency on a scale from 0 to 1 where 0 is completely transparent.

tm_shape(dat_map) + 
  tm_fill("T",
          style = "quantile",
          palette = "Blues"
    ) + 
  tm_borders(alpha=.4) +
  tm_layout(
    legend.outside = TRUE,
    frame = FALSE)

7.12 Add compass

tm_shape(dat_map) + 
  tm_fill("T",
          style = "quantile",
          palette = "viridis"
    ) + 
  tm_borders(alpha=.4) + 
  tm_compass(type = "8star", 
              position = c("RIGHT", "BOTTOM"),
              show.labels = 2,
              text.size = 0.35) +
  tm_layout(legend.outside = TRUE,
            frame = FALSE)

The shape dat_map is invalid. See sf::st_is_valid

Documentation

7.13 Edit the layout

tm_shape(dat_map) + 
  tm_fill("T",
          palette = "viridis",
          style = "quantile",
          title = "Terrorist Attacks"
          ) + 
tm_borders(alpha=.4) + 
tm_layout(
  legend.text.size = 0.7,
  legend.title.size = 1,
  legend.position = c("right", "bottom"),
  legend.outside = TRUE,
  frame = FALSE
  ) 

Documentation

7.14 Interactive map

tmap_mode("view")

tmap mode set to interactive viewing

tm_shape(dat_map) + 
  tm_fill("T",
          palette = "viridis",
          id="province",
          popup.vars=c("T", "growth", "base", "E", "G")
          ) + 
    tm_legend(outside=TRUE) +
  tm_layout(frame = FALSE) 

The shape dat_map is invalid. See sf::st_is_valid

tmap_mode("plot")

tmap mode set to plotting

Documentation

8 Save a new shapefile

st_write(dat_map, "dat_map.shp")

9 References

• https://rdrr.io/cran/tmap/

• Tennekes, M., 2018, tmap: Thematic Maps in R, Journal of Statistical Software, 84(6), 1-39

• https://marcinstepniak.eu/post/interactive-choropleth-maps-with-r-and-tmap-part-i/

• Practical 5: Making maps in R

• Modern Geospatial Data Analysis with R

• https://tlorusso.github.io/geodata_workshop/tmap_package

• https://rstudio-pubs-static.s3.amazonaws.com/229482_b5b182f650434a36ace4eeef2d2bcd83.html#17

• https://mran.microsoft.com/snapshot/2017-12-11/web/packages/tmap/vignettes/tmap-nutshell.html

• Chapter 8 Solutions

• Geocomputation with R website

• https://bookdown.org/lexcomber/brunsdoncomber2e/

10 Datasets

• Non-spatial
• Spatial
• Combined

END